# KMO Kaiser-Meyer-Olkin Measure of Sampling Adequacy
# Function by G. Jay Kerns, Ph.D., Youngstown State University (http://tolstoy.newcastle.edu.au/R/e2/help/07/08/22816.html)
kmo = function( data ){
library(MASS)
X <- cor(as.matrix(data))
iX <- ginv(X)
S2 <- diag(diag((iX^-1)))
AIS <- S2%*%iX%*%S2                      # anti-image covariance matrix
IS <- X+AIS-2*S2                         # image covariance matrix
Dai <- sqrt(diag(diag(AIS)))
IR <- ginv(Dai)%*%IS%*%ginv(Dai)         # image correlation matrix
AIR <- ginv(Dai)%*%AIS%*%ginv(Dai)       # anti-image correlation matrix
a <- apply((AIR - diag(diag(AIR)))^2, 2, sum)
AA <- sum(a)
b <- apply((X - diag(nrow(X)))^2, 2, sum)
BB <- sum(b)
MSA <- b/(b+a)                        # indiv. measures of sampling adequacy
AIR <- AIR-diag(nrow(AIR))+diag(MSA)  # Examine the anti-image of the correlation matrix. That is the  negative of the partial correlations, partialling out all other variables.
kmo <- BB/(AA+BB)                     # overall KMO statistic
# Reporting the conclusion
if (kmo >= 0.00 && kmo < 0.50){test <- 'The KMO test yields a degree of common variance unacceptable for FA.'}
else if (kmo >= 0.50 && kmo < 0.60){test <- 'The KMO test yields a degree of common variance miserable.'}
else if (kmo >= 0.60 && kmo < 0.70){test <- 'The KMO test yields a degree of common variance mediocre.'}
else if (kmo >= 0.70 && kmo < 0.80){test <- 'The KMO test yields a degree of common variance middling.' }
else if (kmo >= 0.80 && kmo < 0.90){test <- 'The KMO test yields a degree of common variance meritorious.' }
else { test <- 'The KMO test yields a degree of common variance marvelous.' }
ans <- list( overall = kmo,
report = test,
individual = MSA,
AIS = AIS,
AIR = AIR )
return(ans)
}
library(psych)
library(GPArotation)
# Load data
raqData <- read.delim(file.choose(), header = TRUE)
# Initial analysis
raqMatrix <- cor(raqData)
cortest.bartlett(raqData)
kmo(raqData)
# PCA with all components
pc1 <- principal(raqData, nfactors = 23, rotate = "none")
pc1
# Scree plot
plot(pc1$values, type = "b")
# PCA with reduced components
k <- 4
pc2 <- principal(raqData, nfactors = k, rotate = "none")
pc2
# Assessing residuals
residuals <- factor.residuals(raqMatrix, pc2$loadings)
residuals
residuals.offdiag <- as.matrix(residuals[upper.tri(residuals)])
hist(residuals.offdiag)
# Orthogonal rotation
pc3 <- principal(raqData, nfactors = k, rotate = "varimax")
pc3
print.psych(pc3, cut = 0.3, sort = TRUE)
# Oblique rotation
pc4 <- principal(raqData, nfactors = k, rotate = "oblimin")
pc4
print.psych(pc4, cut = 0.3, sort = TRUE)
